Density functional theory investigation of the contributions of π-π stacking and hydrogen bonding with water to the supramolecular aggregation interactions of model asphaltene heterocyclic compounds

Context A complex supramolecular process involving electrostatic and dispersion interactions and asphaltene aggregation is associated with detrimental petroleum deposition and scaling that pose challenges to petroleum recovery, transportation, and upgrading. The homodimers of seven heterocyclic model compounds, representative of moieties commonly found in asphaltene structures, were studied: pyridine, thiophene, furan, isoquinoline, pyrazine, thiazole, and 1,3-oxazole. The contributions of hydrogen bonding involving water bridges spanning between dimers and π-π stacking to the total interaction energy were calculated and analyzed. The distance between the planes of the aromatic rings is correlated with the π-π stacking interaction strength. All the dimerization reactions were exothermic, although not spontaneous. This was mostly modulated by the strength of the hydrogen bond of the water bridge and the π-π stacking interaction. Dimers bridged by two water molecules were more stable than those with additional water molecules or without any water molecule in the bridge. Energy decomposition analysis showed that the electrostatic and polarization components were the main stabilizing terms for the hydrogen bond interaction in the bridge, contributing at least 80% of the interaction energy in all dimers. The non-covalent interaction analysis confirmed the molecular sites that had the strongest (hydrogen bond) and weak (π-π stacking) attractive interactions. They were concentrated in the water bridge and in the plane between the aromatic rings, respectively. Methods The density functional ωB97X-D with a dispersion correction and the Def2-SVP basis set were employed to investigate supramolecular aggregates incorporating heterocycles dimers with 0, 1, 2, and 3 water molecules forming a stabilizing bridge connecting the monomers. The non-covalent interactions were analyzed using the NCIplot software and plotted as isosurface maps using Visual Molecular Dynamics. Supplementary Information The online version contains supplementary material available at 10.1007/s00894-024-05922-3.


Introduction
In the last decades, the efficient extraction and transport, marine and terrestrial, of crude oil have become an increasing challenge for the petroleum industry, mainly due to organic and inorganic depositions [1,2].The formation of insoluble aggregates occurs on the surfaces of practically all oil production and transportation systems: in pipelines, reservoirs, containment screens, submerged or surface installations that significantly reduce the recovery of oil [3,4].This generates a large financial loss for the oil industry that, as preventive or remedial measures, must stop production for cleaning or maintenance to increase the useful life of the equipment and the continuity of oil recovery [4][5][6].
Asphaltenes, the densest and most polar fraction of petroleum, are considered the main generators of organic deposition in oil production and transportation systems [7,8].They are a mixture of large molecular weight organic molecules containing polyaromatic chains, side alkyl groups, polar functional groups, and heterocycles with O, N, and S as heteroatoms [9,10].Regarding the molecular structure, asphaltenes can be considered belonging to the continental model, containing a central aromatic moiety with more than seven fused rings and pendant aliphatic chains, or the archipelago model, formed by smaller aromatic moieties interconnected by bridges of alkyl groups [10,11].Experimental and computational studies show that these structures exhibit a strong aggregation preference in solution, forming insoluble solids that reduce oil recovery [12][13][14][15][16][17][18][19][20].Asphaltene aggregation is also detrimental in the context of spill response in the event of accidental oil spills in water [21,22].
Due to the complex constitution of asphaltenes, the insoluble aggregates are stabilized by dispersive, polarization, and electrostatic interactions, arising mainly from London forces, π-π stacking, exchange-repulsion contributions, hydrogen bonding, and dipole-dipole interactions [23][24][25][26].Recently, Hassanzadeh and Abdouss [27], based on studies using the supramolecular assembly model by Gray et al. [28] and the nanoaggregate model by Yen-Mullins [29], have proposed a supramolecular organization model of asphaltene aggregation that combines cooperative binding by acid-base interactions, hydrogen bonding, π-π stacking, and metal coordination, as well as the formation of hydrophobic pockets, porous networks, and host-guest complexes.All these diverse supramolecular interactions are important for the attraction between asphaltene molecules that leads up to generate organic scaling.Based on the supramolecular assembly model [28], Gray et al. have developed organic molecules that would experimentally simulate the aggregation behavior of asphaltenes [30].Figure 1 shows the structure of one of these synthetic asphaltene molecules that contain a bipyridine tethered with ethyl groups to two pyrene moieties.
In previous works [31,32], we investigated computationally the interaction between model asphaltene compounds, reported by Gray et al., to rationalize the experimental information that water traces could enhance the aggregation behavior of asphaltenes.We showed that water molecules could interact with the N atoms of the stacked dimers of asphaltene model compounds, such as PBP and its analogs with varied aromatic hydrocarbon groups.The interaction strength comparison between the π-π stacking and hydrogen bond showed that these interactions have almost the same contribution to the stabilization of the nanoaggregate [31,32].
In the present study, we explore the formation of bridges of water molecules spanning between the stacked homodimers of N-, O-, and S-containing heterocycles and evaluate the enhancement of the supramolecular interaction strength.The heterocycles pyridine, thiophene, isoquinoline, and furan were selected due to their presence as moieties in asphaltene structures [33][34][35][36][37][38].Additionally, the heterocycles pyrazine, 1,3-oxazole, and 1,3-thiazole were included to study aromatic compounds with two heteroatoms.The structures of these compounds are shown in Fig. 2. The contributions of π-π stacking and hydrogen bonding interaction to the aggregate stabilization were investigated using non-covalent interaction and energy decomposition analysis (EDA).The interaction energy was rationalized in terms of geometric and energetic parameters inherent to the dimers and monomers.

Computational methods
The DFT calculations were performed with the ωB97X-D exchange-correlation density functional with a dispersion energy correction [39].Several studies have previously shown successful results with the account of the interaction of π-conjugated oligomers [40], aggregation of asphaltene model compounds [19,32], and the interaction of polycondensed aromatic molecules [41].We employed the Def2-SVP basis set for all atoms.This is a split-valence double-zeta basis set with polarization functions for all atoms, proposed by Ahlrichs and Weigeng [42].All calculations were performed in Gaussian09 [43].After full geometry optimization, the second-order force constant matrix was calculated to confirm that the optimized geometry is a genuine minimum on the potential energy surface.The thermodynamic results are important for the computational determination of the enthalpy and Gibbs free energy (at 1 atm and 298 K) of the formation of supramolecular aggregates.The basis set superposition error (BSSE) was accounted for by the counterpoise correction procedure [44].The B3LYP/Def2-SVP method was employed for the energy decomposition analysis using the GAMESS software [45,46].The EDA procedure [47][48][49] decomposes the total interaction energy into five components: electrostatic (E Elec ), polarization (E Pol ), exchange (E Xc ), dispersion (E Disp ), and Pauli repulsion (E Pauli ).The sum of the polarization and exchange terms yields the covalent component of the interaction.
Non-covalent interactions between the monomers in the dimer structure and between the monomers and water molecules in the aggregate structure were also analyzed.For this, we used the NCIplot software [50].This software calculates the attractive and repulsive interactions present in the system as a function of the electron density and its reduced gradient.Non-covalent interactions are plotted as an isosurface map using Visual Molecular Dynamics [51].

Geometry optimization
The geometries of 28 supramolecular aggregates containing homodimers and water molecules were fully optimized with the ωB97X-D/Def2-SVP method.The optimized structures for the aggregates containing dimers with and without water molecules are shown in Fig. 3.For pyridine, furan, thiophene, and isoquinoline which have one heteroatom, only one water bridge is formed, whereas for the pyrazine, oxazole, and thiazole rings that have two heteroatoms, two water bridges are formed, stabilizing the dimer structure.
For pyridine, furan, thiophene, isoquinoline, oxazole, and thiazole dimers, the bridge can be formed by one, two, or three water molecules, with the aromatic moieties retaining almost parallel orientation.For the pyrazine dimer, only bridges with one or two water molecules are formed.The pyrazine aggregate, composed of three water molecules (Pyra 2 (H 2 O) 6 ), upon optimization converges to a structure where only two water molecules are in the bridge with the third water molecule being hydrogen bonded to the water bridge.As pyrazine has two heteroatoms, two water bridges may be formed, making the aggregate more rigid and keeping the interplanar distance short because of the stabilizing effects of the water bridges.The high symmetry of the pyrazine ring strengthens the π-stacking interaction and does not allow the aromatic rings to be sufficiently far from each other to accommodate a third water molecule in the bridge.
For isoquinoline, both the dimer without any water molecule (Iso 2 ) and the aggregate with three water molecules (Iso 2 (H 2 O) 6 ) have aromatic rings that are in nonparallel planes, whereas the aggregates of isoquinoline with one (Iso 2 (H 2 O) 2 ) or two (Iso 2 (H 2 O) 4 ) water molecules in the bridge have well-aligned aromatic planes.The bridge of one or two water molecules forces the aromatic planes to be one above the other.The dimer without water (Iso 2 ) does not exhibit this effect.The aggregate with three water molecules has more degrees of freedom and a longer bridge that causes the distortion of the heterocyclic plane rather than constraining those to a parallel configuration.
The thiazole aggregates with one or two water molecules in the bridge exhibit more substantial horizontal displacement than the structures with three water molecules (Thia 2 (H 2 O) 6 ) and the water-free dimer (Thia 2 ).This is related to the interaction of the water molecules with the heteroatom of the aromatic system and the resultant competition between π-π stacking and electrostatic interaction.Moreover, in the aggregates with one (Thia 2 (H 2 O) 2 ) or two (Thia 2 (H 2 O) 4 ) water molecules per bridge as well as in the water-free dimer (Thia 2 ), the heterocycles are in a staggered configuration; i.e., the N and S atoms of one of the monomers are on opposite sides from those in the other monomer.Thus, the water bridges span between different heteroatoms, forming N⋯H 2 O⋯S networks, whereas for the bridge containing three water molecules, the interaction is with the same heteroatom, forming N⋯H 2 O⋯N and S⋯H 2 O⋯S hydrogen bonding networks.We also calculated the thiazole aggregates with bridges containing one (Thia 2 (H 2 O) 2 ) and two (Thia 2 (H 2 O) 4 ) water molecules involving the same heteroatom, i.e., N⋯H 2 O⋯N and S⋯H 2 O⋯S, but these are less stable than the Thia 2 (H 2 O) 2 and Thia 2 (H 2 O) 4 shown in Fig. 3.As S has a larger atomic radius than N, its orbitals do not align with the one of the other atoms of the aromatic ring and the π-stacking interaction is more effective by the side of the N atom.This is the reason why the aromatic rings are horizontally displaced, as in the dimer without water molecules (Thia 2 ) the S atoms are in positions opposite to each other.In the aggregate containing three water molecules per bridge (Thia 2 (H 2 O) 6 ), the longer bridge allows the water trimer to interact with the same heteroatoms.It was not possible to optimize the aggregate with a water trimer bridge spanning between different heteroatoms.

Energy calculation
The stabilization of the dimer structures is due to the π-π stacking and the water bridge hydrogen bond interactions.The π-π stacking interaction between the rings in each dimer structure was accounted for by means of the ASM method (activation strain model).This model, proposed by Fernández and Bickelhaupt [52], is a computational approach that uses electronic structure calculations to rationalize the factors that control stability in each stationary point along a reaction coordinate.According to the ASM model, the binding energy (E BIND ) can be decomposed into two contributions along the reaction coordinate (Eq.1).The first one is the strain or distortion energy, ΔE DEF , related to the deformation of the reactants, mainly affected by the rigidity of their molecular structure and by distortion in pendant groups.The second component is the interaction energy, ΔE INT , between the reactants, related to the binding capacity between the reactants, which, in our case, accounts to the π-π stacking interaction.
In general, ΔE DEF is a positive term, which destabilizes the system, and the ΔE INT is a negative term, stabilizing the system.The ΔE INT can be decomposed into stabilizing electrostatic interaction (coulomb) between fragments ∆V elst ; destabilizing interaction derived from the overlap of filled orbitals ∆E Pauli ; orbital interaction energy ∆E oi , responsible for the charge transfer (HOMO-LUMO interactions, for example) and polarization (mixing unoccupied and occupied orbitals in the different fragments); and the interaction due to dispersion forces ∆E disp .These data are presented in the Supplementary Material.Table 1 presents the values of ΔE INT between the heteroaromatic ring dimers and the interplanar distance of the structures shown in Fig. 3.
The analysis of Table 1 shows that the strongest ΔE INT interaction, corresponding to the lowest energy value, is for the aggregates without any water molecule in the structure.The larger freedom or reduced constraints of the π-electronic clouds of the rings that can undergo horizontal displacement enhance π-π stacking interaction.The largest value of ΔE INT is for the isoquinoline ring, almost twice the values found for the other rings; it is followed by the pyrazole and pyridine rings.The isoquinoline is the only system with two fused rings, which enhances π-π stacking interaction.Thus, for water-free dimers, ΔE INT is modulated by the size of the aromatic system.
The structures linked by water molecules are more rigid and constrained, which cannot assume the alignment between the ring planes required to accommodate the electron density for the π-stacking interaction.In all cases, we note that the aggregates containing a bridge of just one water molecule have the smallest (less negative) ΔE INT interaction among the water-bridged dimers.The bridge formed by only one water molecule makes the system more rigid, restricting the planes from the adequate horizontal displacement for optimal π-stacking interaction.Generally, the aggregates with two water molecules in the bridge have the second lowest ΔE INT , followed by the dimers with three water molecules per bridge.The bridges containing two water molecules apparently provide more favorable structural arrangements for the aromatic planes that enhance the π-stacking interaction than the bridges with three water molecules.
In Table 1, we also show the interplanar distance in the dimers.For the pyridine and furan dimers, we notice that as the number of water molecules in the bridge increases, ( 1) the distance becomes larger.This is due to the hydrogen bond interaction between the dimers and the water bridge that is so effective that enhances the distance between the aromatic planes.For the thiophene and pyrazole aggregates with bridges of two water molecules (Thio 2 (H 2 O) 2 and Pyra 2 (H 2 O) 4 ), the interplanar distance is smaller than that for the aggregates with bridges of three water molecules.This is attributable to the more effective accommodation of the aromatic electronic cloud overlap between the planes.For the oxazole and isoquinoline aggregates, the dimers with two (Oxa 2 (H 2 O) 2 and Iso 2 (H 2 O) 2 ) and three (Oxa 2 (H 2 O) 3 and Iso 2 (H 2 O) 3 ) water molecules in the bridge have shorter interplanar distances than those with one or without any water molecule in the bridge.For the isoquinoline system, the aromatic planes deviate from parallel to a crossed configuration as the number of water molecules increases, disrupting the π-π stacking interaction; however, the shorter interplanar distance likely partially offsets for the decreased parallel alignment.The thiazole dimer with two water molecules (Thio 2 (H 2 O) 4 ) in the bridge has the longest interplanar distance of all aggregates analyzed.The π-π stacking interaction and dipole-dipole interaction of the S atom with the water molecules in the bridge are not as strong as the hydrogen bonding within the bridge.
The enthalpy (ΔH) and Gibbs free energy (ΔG) for the supramolecular aggregation were calculated considering either none or a cluster of 1-3 water molecules (Eqs. 2 and 3).
where H(orG) complex is the enthalpy (or Gibbs free energy) of the aggregate with or without water molecules in the bridge, H(orG) monomer is the enthalpy (or Gibbs free energy) of each monomer, and E (H 2 O) x is the enthalpy (or Gibbs free energy) of one ( x = 1 ), cluster of two ( x = 2 ), or cluster of three ( x = 3 ) water molecules.
In Table 1, we show the π-π stacking interaction and thermodynamic results for all the structures shown in Fig. 3.In the water-free dimers, the ΔH and ΔG are stabilized only due to the π-π stacking interaction between the monomers.In aggregates with the water molecules, in addition to the π-π stacking interaction, the hydrogen bonds also help stabilize the complexes.
The analysis of Table 1 shows that the aggregations are exothermic (negative ΔH), although not spontaneous (positive ΔG).In general, the aggregates containing two or three water molecules per bridge are the ones that have the smallest values of ΔH, suggesting that the complexes formed with two water molecules are the most stable, followed by those with three water molecules.This result has previously been reported based on the calculation of the aggregation of model asphaltene compounds [32].As the ΔH analysis accounts for both the π-stacking and the hydrogen bond terms, the aggregates with two water molecules per bridge must have the strongest stabilization by hydrogen bonding.The ΔG analysis shows that the structures without any water molecules have the smallest values, although still being positive, in agreement with the experimental report that these molecules do not aggregate spontaneously [31].
The water bridges increase the distance between the heterocycles, by a maximum of 0.1 Å, weakening the π-π stacking interaction; however, the stabilization of the aggregate seems to be compensated by the hydrogen bonds, since the formation of complexes with two or three water molecules is exothermic.

Intermolecular interactions
To examine the non-covalent interactions between the monomers in the dimer structure and between the monomers and water bridges, the NCI plots of the optimized structures of the supramolecular aggregates are presented in Fig. 4. The isosurface colors vary between green and blue colors.The green color represents a weak favorable non-covalent interaction, such as the van der Waals interaction.The blue color represents strong favorable noncovalent interactions, such as conventional hydrogen bonds.Unfavorable and repulsive interactions, represented in red color, are not observed on the isosurfaces of any of the aggregates investigated.
The analysis of Fig. 4 shows that the green area corresponds to the π-π stacking interaction between the heteroaromatic rings.The green color indicates that this interaction is weak, as expected for this type of interaction, and the dispersed isosurface shows a delocalized interaction.There is no significant difference in the plotted area for the aggregates with one, two, or three waters in relation to the dimers without water.As isoquinoline is composed by two fused aromatic rings, its π-stacking interaction is more diffuse and is seen with an enlarged green area.Light blue areas, representing strong attractive interactions due to hydrogen bonding, are also seen in Fig. 4. The blue regions represent localized interactions and are not as dispersed as green ones.Furthermore, as the number of water molecules per bridge increases, the localized interactions shown in blue, corresponding to hydrogen bonds, increase, indicating a stronger bridging interaction, as expected for hydrogen bonding networks.For the interaction between the S atom and the water molecules, in thiophene and thiazole, we could not observe blue areas, indicating weaker dipole-dipole interactions.
The analysis of Table 2 shows a trend that is observed for all the components of the interaction as well as for the total interaction energy.As a general trend, the total interaction energy increases when increasing the number of water molecules in the bridge.However, the incremental difference is much more relevant for the first water molecule than for the second or the third one.For example, the difference in the E tot values for the aggregates with one water molecule per bridge to those with two water molecules per bridge is 9.15 ± 3.55 kcal mol −1 , whereas the difference for the aggregates with two water molecules per bridge to those with three water molecules per bridge is only 1.05 ± 1.65 kcal mol −1 .This shows that the aggregates with two and three water molecules are significantly more stabilized than those with only one water molecule per bridge.Also, the energies of aggregates with two or three water molecules in the bridge do not vary substantially.
The dispersion term (Table 2), accounting for long-range interactions, has the smallest variation (standard deviation of 3.68 kcal mol −1 ).It also changes more strongly from the aggregates with one water molecule per bridge to the ones with two water molecules per bridge than for additional water molecules.The structures with the stronger electrostatic term also have the largest repulsion term (E Pauli ).
In Fig. 5, we present the electrostatic and covalent components of the total interaction energy.In this model, the ionic character of the interaction is accounted for by the E Elec term, which comes mainly from opposite charge attraction sites.The covalent component is due to the sum of the E Pol and E Xc terms and considers the overlap of the atomic orbitals that compose the interaction.We can see for all dimers that the covalent character is almost two times larger than the ionic character, corresponding to a stabilization of − 12.54 ± 4.41 kcal mol −1 .As the pyrazole, thiazole, and pyrazine aggregates have two water bridges, the stabilization per water bridge is the total value divided by 2. Considering the stabilization energy per water bridge, the most stable aggregates are those with hydrogen bonds between the heteroatom of the dimer and water of the bridge, i.e., interaction of H of water bridge with the O or N atom of the heterocycle.The S-containing heterocycles thiophene and thiazole with dipole-dipole interaction between the H atom of the water bridge and S atom of the heterocycle have lower stabilization energy.

Conclusion
We investigated the interaction of the supramolecular aggregation of seven heterocyclic aromatic compounds as waterfree dimers as well as dimers with water bridges bonded to the heteroatoms and spanning between the organic planes.We observed, for most of the aggregates, that the interaction is favored by bridges composed of two water molecules.Only for the 1,3-thiazole the favorite bridge has three water molecules, probably due to the softness of the sulfur atom.The π-π stacking interaction analysis showed that the waterfree dimers have the strongest interaction, followed by the dimers with two water molecules per bridge.In almost all cases, the π-π stacking interaction strength is modulated by the interplanar distance between the monomers in the dimer structure.The ΔH analysis showed that aggregation is an exothermic process.The most stable aggregate for each heterocycle is the system with two water molecules per bridge.The ΔG analysis showed nonspontaneous aggregation processes with the smallest values for the dimer without any  water molecule in the bridge.The NCI plot analysis identified strong interaction sites around the water molecules, representing hydrogen bonding interactions, and weak attraction between the planes of the organic molecules, representing the π-stacking interactions.The hydrogen bonds between the dimers and the water bridges were decomposed by using the EDA method.The results indicated that the covalent character (polarization and exchange) of the interaction is almost twice as large as the electrostatic term.We also noticed that one water molecule in the bridge led to a small stabilization of the aggregate, whereas two or three water molecules in the bridge add a considerable stabilization to the supramolecular system, with the aggregates having two water molecules per bridge being the most stable.Our findings justify the conclusion that the π-π stacking interaction is as important as hydrogen bonding for the stabilization of the dimers bridged by water molecules.

Fig. 4
Fig. 4 Optimized structures of the dimers without water, and with one, two, or three water molecules per bridge, showing the non-covalent interaction obtained with NCIplot

Table 1 ΔE
INT(in kcal mol −1 ), ΔH (in kcal mol −1 ), and ΔG (in kcal mol −1 ) for the formation of the aggregates presented in Fig.3.The distance (D) between the centers of the aromatic rings is in angstrom (Å).Values for ΔE INT , ΔH, and ΔG are corrected by BSSE *ΔE INT , ΔH, and ΔG values were calculated taking the isolated monomers and one, a cluster of two, or a cluster of three water molecules as reference **The distance shown is the distance between the centers of the aromatic rings

Table 2
The EDA components E tot , E Elec , E Pol , E Xc , E Disp , and E Pauli in kcal mol−1